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A quantum-classical limit of the canonical equilibrium time correlation function for a quantum 
system is derived. The quantum-classical limit for the dynamics is obtained for quantum systems 
comprising a subsystem of light particles in a bath of heavy quantum particles. In this limit the 
time evolution of operators is determined by a quantum-classical Liouville operator but the full equi- 
librium canonical statistical description of the initial condition is retained. The quantum-classical 
correlation function expressions derived here provide a way to simulate the transport properties 
of quantum systems using quantum-classical surface-hopping dynamics combined with sampling 
schemes for the quantum equilibrium structure of both the subsystem of interest and its environ- 
ment. 



I. INTRODUCTION 

The dynamical properties of systems close to equilibrium may be described in terms of equilibrium time correlation 
functions of dynamical variables or operators. For a quantum system with Hamiltonian H at temperature T with 
volume V, linear response theory shows that the time correlation function of two operators A and B, needed to obtain 
transport properties, has the Kubo transformed form 
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where (3 — (fe^T) -1 , p e = Z~^e~@ H is the quantum canonical equilibrium density operator, Zq = Tre _/9ff is the 
canonical partition function and, in the second line, t\ = t — ih(f3 — A) and £2 = t — ih\. The evolution of the operator 
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A(t) is given by the solution of the Heisenberg equation of motion, dA(t)/dt — i[H, A(t)], where the square brackets 
denote the commutator. 

While such correlation functions provide information on the transport properties of the system, their direct compu- 
tation for condensed phase systems is not feasible due to our inability to simulate the quantum mechanical evolution 
equations for systems with a large number of degrees of freedom. While approximate schemes have been devised 
to treat quantum many body dynamics, for example, quantum mode coupling methods have proved useful in the 
calculation of collective modes for some applications [3j, we are primarily concerned with methods that approximate 
the full many body evolution of the microscopic degrees of freedom. In many circumstances only a few degrees of 
freedom need to be treated quantum mechanically (quantum subsystem) while the remainder of the system with which 
they interact can be treated classically (classical bath) to a good approximation. Examples where such a description 
is appropriate include proton and electron transfer processes occurring in solvents or other chemical environments 
composed of heavy atoms. Quantum-classical methods have been reviewed by Egorov et al. 4] and one form of a 
quantum-classical approximation has been assessed in the weak coupling limit where there is no feedback between the 
quantum and classical subsystems. Although it is difficult to determine transport properties such as the reaction rate 
constant from the full quantum time correlation function when the entire system is treated quantum mechanically, 
methods are being developed to carry out such calculations. 5] Mixed quantum-classical methods also provide a route 
to carry out nonadiabatic rate calculations. 

DP 

A number of schemes have been proposed for carrying out quantum dynamics in classical environments. |y, |7 

We focus on approaches where the evolution is described by a quantum-classical Liouville equation 
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For a quantum system coupled to a classical environment it is possible to derive an 



evolution equation for dynamical variables or operators (or the density matrix) by an expansion in a small parameter 
that characterizes the mass ratio of the light and heavy particles in the system. The quantum-classical analog of the 
Heisenberg equation of motion is, 



-A w (X,t) = ±[6 w ,Aw(t)]-^[\lI\\.Au(h\-{A n [t).II [V } 



= iCA w (X,t) 
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Here Aw(X, t) is the partial Wigner representation |2l|,|22j of a quantum operator; it is still an operator in the Hilbert 
space of the quantum subsystem but a function of the phase space coordinates X = (R, P) of the classical bath. In 
this equation {•••,•••} is the Poisson bracket and C is the quantum-classical Liouville operator. A few features 
of quantum-classical Liouville dynamics are worth noting. This equation of motion includes feedback between the 



classical and quantum degrees of freedom. The environmental dynamics is fully classical only in the absence of 
coupling to the quantum subsystem. In the presence of coupling the environmental evolution cannot be described 
by Newtonian dynamics, although the simulation of the quantum-classical evolution can be formulated in terms of 
classical trajectory segments. |2l| For harmonic environmental potentials with bilinear coupling to the quantum 
subsystem the evolution is equivalent to the fully quantum mechanical evolution of the entire system. Quantum- 

n 

classical simulations of the spin-boson model are in accord with the numerically exact quantum results 23] and have 
been used to test quantum-classical simulation algorithms. [24I 2{j | 

Equation (J2J is valid in any basis and an especially convenient basis for simulating the evolution by surface-hopping 
schemes is the adiabatic basis, {\ai; R >}, the set of eigenstates of the quantum subsystem Hamiltonian in the presence 
of fixed classical particles. In this case the matrix elements of an operator A^" 1 (X,t) =< a\\ R\A W (X, t)\a[; R > 
satisfy 

±A^(X,t) = i J2 C aia[Pl0i A^(X,t) , (3) 

0i0[ 
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where C aia > p x p> denotes the representation of the quantum-classical Liouville operator in the adiabatic basis. 
This equation may be solved using surface-hopping schemes that combine a p robabilistic description of the quantum 
transitions interspersed with classical evolution trajectory segments. [2^. LJ, I2II I27 . 2$j Although further algorithm 
development is needed to carry the simulations to arbitrarily long times, the quantum-classical evolution is not a 
short time approximation to full quantum evolution. Rather, it is an approximation to the full quantum evolution 
for arbitrary times since the quantum-classical evolution is derived at the level of the Liouville operator and not the 
quantum propagator. Given the evolution equation J2J) (and the corresponding quantum-classical Liouville equation 
for the density matrix) one may construct a statistical mechanics of quantum-classical systems 2, 3^| 
transport properties such as chemical rate constants i^] 



and compute 



from the correlation functions obtained from this analysis. 
In this article we consider another route to determine quantum-classical correlation functions for transport proper- 
ties. We begin with the full quantum mechanical expression for the time correlation function (Eq. ijTJl) and take the 
limit where the dynamics is determined by quantum-classical evolution equations for the spectral density that enters 
the correlation function expression. While the calculations leading to our expression for the correlation function are 
somewhat lengthy, the final result has a simple structure: 
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This expression for the time correlation function retains the full quantum statistical character of the initial condition 
through the spectral density function W (Eq. H4()|l below) but the forward and backward time evolution of the 
operators b]^ 1 ^ 1 and A^ 2 , respectively, is given by the solutions of the quantum-classical evolution equation 
Consequently, one may combine algorithms for determining quantum equilibrium properties with surface-hopping 
algorithms for quantum-classical evolution to estimate the value of the correlation function. Quantum effects enter in 
all orders in this expression for the correlation function. In addition to the fact that the initial value of the spectral 
density contains the full quantum equilibrium statistics, since the quantum-classical Liouville operator appears in 
the exponent in the propagator, the quantum-classical propagator contains all orders of h, albeit in an approximate 
fashion. 

The outline of the paper is as follows: In Sec.[n]we construct the partial Wigner representation of the quantum time 
correlation function and obtain expressions for the spectral density and its matrix elements in an adiabatic basis. In 
Sec. IIIII we derive a quantum-classical evolution equation for the matrix elements of the spectral density and establish 
a connection to quantum-classical Liouville evolution. The results of these two sections are used in Sec. UVI to obtain 
Eq. J3J. In Sec. we analyze the initial value of the spectral density in the high temperature limit. The conclusions 
of the study are given in Sec. IVII while additional details of the calculations are presented in the Appendices. 

II. PARTIAL WIGNER REPRESENTATION OF QUANTUM CORRELATION FUNCTION 

We consider quantum systems whose degrees of freedom can be partitioned into two subsets corresponding to light 
(mass m) and heavy (mass M) particles, respectively. We use small and capital letters to denote operators for phase 
points in the light and heavy mass subsystems, respectively. In this notation the Hamiltonian operator for the entire 
system is the sum of the kinetic energy operators or the two subsystems and the potential energy of the entire system, 
H = P 2 /2M + p 2 /2m + V(q, Q). 

We are interested in the limit where the dynamics of the heavy particle subsystem is treated classically and the light 
particle subsystem retains its full quantum character. To this end it is convenient to take a partial Wigner transform 
of the heavy degrees of freedom and represent the light degrees of freedom in some suitable basis. 

In order to carry out this program we begin with the quantum mechanical Kubo transformed correlation function 
and write the trace over the heavy subsystem degrees of freedom in the second line of Eq. using a {Q} coordinate 



representation, 
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C AB (t;(3) = 4/ dX-^-Tr' ifldQi < Qi|£t|Q 2 >< Q 2 \e* t * 6 \Q 3 >< Q 3 \A\Q 4 >< Q 4 M t2 ^|Qi > . (5) 

The prime in Tr' refers to the fact that the trace is now only over the light particle subsystem degrees of freedom. 
Making use of the change of variables, Q\ = R\ — Zi/2, Q 2 = R\ + Zi/2, Q3 = R2 — Z2/2 and Q4 = R2 + Z2/2, this 
equation may be written in the equivalent form, 

C AB (t;(3) = 1 f 'dX^Tr' [ TidR^ < R, - ^\&\R X + ^-><R 1 + ^-\e^"\R 2 - ^ > 
P Jo Z Q J fj[ 2 1 2 2 
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x < R 2 --±\A\R2 + ^-><R2 + Y-\e-^ H \R l ~^-> . (6) 

The next step in the calculation is to replace the coordinate space matrix elements of the operators with their 
representation in terms of Wigner transformed quantities. The partial Wigner transform of an operator is defined by 

m, 



A W (R 2 ,P 2 )= J dZ 2 e* p * z *(R 2 - ^-\A\R 2 + y> - 



(7) 



while the inverse transform is 



(R 2 y |i|#2 + f ) = J dP 2 e-* p * z *A w (R 2 ,P,) 



(8) 



Here is the dimension of the heavy mass subsystem. The partially Wigner transformed operator A^iX^) is 
a function of the phase space coordinates X2 = (R 2 ,P 2 ) and an operator in the Hilbert space of the quantum 
subsystem. It is convenient to consider a representation of such operators in basis of eigenfunctions. In this paper 
we use an adiabatic basis since, through this representation, we can make connection with surface-hopping dynamics. 
The partial Wigner transform of the Hamiltonian H is H w = P 2 /2M + p 2 /2m + V w (q,R) = P 2 /2M + h w (R). 
The last equality defines the Hamiltonian hw(R) for the light mass subsystem in the presence of fixed particles 
of the heavy mass subsystem. The adiabatic basis is determined from the solutions of the eigenvalue problem, 
hw(R)\at; R >= E a (R)\a; R >. The adiabatic representation of A\y{X2) is, 

A W (X 2 ) = l«2;#2 > A^(X 2 ) < a' 2 ;R 2 \ , (9) 

OL20l-' 2 

where A% a '*(X 2 ) =< a 2 ; R 2 \A W (X 2 )\a' 2 ; R 2 >. 

By inserting Eq. © into Eq. (JSJ we can express the coordinate representation of the operator A as 

< R 2 - ^-\A\R 2 + y>= £ / dP2 z~ iP2 Z2 \^R2 > A^ a \x2) < a' 2 ;R 2 \ . (10) 
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Then, substituting the result into Eq. © (along with a similar representation of the W operator), we obtain, 

C AB (t;f3) = i ^dX J2 J fidX.B^iXMT' 2 ^) 

x W a 'i aia 2 a2 (X 1 ,X 2 ,t;fJ,\) . (11) 

Here we defined the matrix elements of the spectral density by 

Zq J f\ 2 2 

x < a' 2 ;R 2 \ < R 2 + ^e^"^ - ?± > \a 1 ;R 1 > * . (12) 

Our task is now to find an evolution equation for W aiai0l2a2 (Xi, X 2 ,t; f3, A) in the mixed quantum-classical limit. 
Before doing this we observe that the expression for the quantum correlation function in the partial Wigner represen- 
tation is equivalent to an expression involving full Wigner transforms of the operators. 



A. Relation to full Wigner representation 

Since the correlation function is independent of the representation we choose for the light and heavy mass sub- 
systems, we may also represent the light mass subsystem in terms of a Wigner transform instead of a set of basis 
functions. To establish the connection between these two forms of the correlation function we note that the full 
Wigner transform of the the operator A is given by 

A w (x 2 ,X 2 ) = J dz 2 e^ Z2 <r 2 - ^\A w (X 2 )\r 2 + y > , 

= e^ z ^ a2 (r 2 - ^;R 2 )A a ^(X 2 )cb^(r 2 + |;i? 2 ) , (13) 

a 2 a! 2 

where x 2 — {r 2 ,p 2 ) and 4> a2 (r 2 ; R 2 ) =< r 2 \a 2 ; R 2 >. We have used Eq. @ to write the second line of this equation. 
The inverse of this expression is 
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l vT 2 (* 2 ) = 7^^7 / d V2dir 2 - ^)d(r 2 + ^) e^^ 2 A w (x 2 , X 2 )^(r 2 - ^;i? 2 )^(r 2 + ^;R 2 ) , (14) 



where is the dimension of the light mass subsystem. Inserting this expression for A^ a2 (X 2 ) (and the analogous 
expression for B^ 1 " 1 (Xi)) into Eq. (|ll|l we find 



C AB (t,l3) 



\ [ dX f fldxidXi B j w (x 1 ,X 1 )A w (x 2 ,X 2 ) 
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xW(x 1 ,X 1 ,x 2 ,X 2 ;t,0,X) , (15) 



where 



X^(r 2 + |;i? 2 )C 2 (r 2 - |; iZ 3 )W^°i"»°" X 2 , t; 0, A) * . (16) 

Using the definition of the matrix elements of W in Eq. (|12|) and performing the sums on states we may write this 
equation in the equivalent form, 



W{x 1 ,x 2 ,X 1 ,X 2 ,t;(3,X) = -J- [fXdz^e-^^+^e-i 
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x < ri + z ±\ < Rl + ?±\ e it-H\ R2 _| >h _| > 

x<r2 + | l<fl2 + | |e -*«,* |fll -| >|ri _|>_J_, (17) 

where v = vi + Vh. Equation l|17|) gives the spectral density in the full Wigner representation while Eq. I|16(l relates 
this quantity to its matrix elements in the light mass subsystem basis. In particular, letting W be a (super) matrix 
whose elements are Vl^" 1 " 1 " 2 " 2 , Eq. I|16l) may be written formally as 

W = ToW, (18) 

where T is the transformation specified by Eq. (|lfci|) . 

For future use, we note that the inverse of this expression is given by 

f 2 

W a '^ a '^(X 1 ,X 2 ,t-p,\) = I \[dxdz l e^- z ^ z -U a Ari- Z ^;Ri)^ ri + ^ R ^ 

1—1 " 

x<^ 2 (r 2 - y; Jfc)&/ (r 2 + ^;R 2 )W{xx,x 2 , X u X 2 ,t; /?, A) , (19) 

as can be verified by substituting Eq. (|17fl into Eq. 1|19|) and performing the integrals. Like Eq. i|ltj|) . Eq. I|19|) gives 
a relation between W and its matrix elements but now in the opposite direction, relating VK" 1 " 1 " 2 " 2 to W . Using a 
formal notation like that in Eq. 1)18(1 . we can write Eq. 1)19(1 as 

w = r 1 o w , (20) 

which defines 1 , the inverse of T. 

III. QUANTUM-CLASSICAL EVOLUTION EQUATION FOR SPECTRAL DENSITY 

The quantum-classical evolution equation for matrix elements of the spectral density Vl /a i QlQ 2 Q2 (X\, X 2l t; (3, A) can 
be obtained using various routes. In this paper we first derive a quantum-classical evolution equation for the spectral 
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density W(xi, X2, Xx,X2, t; (3, A) in the full Wigner representation. We then change to an adiabatic basis representation 
of the quantum subsystem to obtain our final result for the evolution equation for W aiaia2a2 {X\, X2, t; (3, A). 

A. Evolution of W in the full Wigner representation 

The evolution equation for X2, Xi, X2, t; f3, A) can be obtained by differentiating its definition in Eq. I|17(l 

with respect to time and then inserting complete sets of coordinate states to obtain a closed equation in W. The 
result was obtained earlier by Filinov et al. 32( and, for our composite system, is given by 

^-W(x 1 ,x 2 ,X 1 ,X2,t; /3, A) = ~{iL^\x u X x ) - ilf\x 2 , X 2 )) W(x 1 ,x 2 ,X 1 ,X 2 , t; f3, A) 

1 f 2 

+ - J W_ds l dS l {uj 1 (r 1 ,s 1 ,Ri,S 1 )8(s 2 )5(S 2 ) - ^2^2, s 2 , R 2 , S 2 )S(si)5(Si)j 

J i=l 

x W{x 1 -ir 1 ,X 1 -n 1 ,X2-n2,X2-n 2 ,t;(3,X) . (21) 
Here we have introduced the notation iTi = (0, Si), IT = (0, Si). The classical free streaming Liouville operators are 

m oqt M oRi 
for (i = 1,2). The Ui functions under the integral are defined by 

uj l (r l ,R i ,s l ,S l ) = 2 I ' dfid^Viri-fi^R^ R\) 
n(iTh) u J 

2 2 

x sin(- Si ■ fi + -Si ■ Ri) . (23) 

B. Scaled equation of motion 

In order to take the quantum-classical limit of Eq. (|21J) . we consider systems for which the ratio between the 
light particle mass and the heavy particle mass is small, and employ the same mass scaling used in Ref. 21 1 to 
obtain the quantum-classical Liouville equation. One is naturally led to consider an expansion in the small parameter 
fi = (m/M) 1 / 2 from the following arguments. Consider a unit of energy eoi say the thermal energy eo = /3 _1 , a unit 
of length A m = h/p m , where p m = (meo) 1 / 2 is the unit of momentum of the light particles, a unit of time to = 
and a unit of momentum for the heavy particles Pm — (M^o) 1 ^ 2 - These units may be used to scale the coordinates 
of the system so that the magnitude of the scaled momentum of the heavy particles, P/Pm, is of the same order of 
magnitude as that for the light particles, p/p m - Only momenta are scaled by different factors; characteristic lengths 
are scaled by the light particle thermal de Broglie wavelength, A m . This is analogous to the scaling used to derive the 
equations of Brownian motion for a heavy particle in a bath of light particles. 
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In the following, we denote scaled quantities with a prime; e.g., r' — r/X m , R' — R/X m , p' — p/pm, P' = P/Pm, 
t' = t/to, etc. The scaled version of the equation of motion for W has the same form as Eq. H21(l but with all quantities 
replaced by their primed dimensionless counterparts. The scaled operators and functions in this equation are 

tip 04Xi)=^A +Mi ^. JL, (24) 

and 

x sin{2s' i -f' i + 2S' i -R' i ). (25) 

In writing the last line of the u>[ equation we have performed the change of variables R\ = in the dummy 

variable in the integration in order to move the /x dependence from the sine factor to the potential, which is more 
convenient for taking the classical limit. We see that the classical free streaming evolution is linear in fi but the 
quantum kernel has all powers of \x. 



C. Quantum-classical equation for W 

For (i << 1 the quantum-classical limit is obtained expanding the evolution operator up to linear terms in /i. Since 
momentum is related to the de Broglie wavelength of a particle, this procedure is equivalent to averaging out the 
short de Broglie oscillations of the heavy particles on the scale of the long de Broglie oscillations of the light particles. 

The expansion of the evolution operator is obtained from the expansion of the interaction potential to linear order 
in the small parameter /j,, 



V'{r[ -f'^Ri- K) - V(r> - f>, R£ - vR> ■ ^ l dR > 1 + ■ 



(26) 



Inserting this expansion in Eq. I|25|) . working out the integrals, substituting the result into the scaled version of Eq. 121|) 
and finally going back to unsealed coordinates (the details are given in Appendix^) one obtains the quantum-classical 
equation for W in unsealed coordinates as 



^-W(x u x 2 ,X 1 ,X 2 ,t;p,X) = -i 



iLf\x 1 )+iL 1 {X l )-iLf\x 2 ) -iL 2 (X 2 ) W{x 1 ,x 2) X 1 ,X 2 ,t;f3,X) 



dsids 2 



1 



dr 



S(s 2 )V(r 1 -f,R x ) sin ( ^-^-) ~ 5( Sl )V(r 2 -f,R 2 ) sin ( 2 " ' ' 



1 

+ 2 



S(s 2 )AF 1 (R 1 , a x ) ■ — - S(s 1 )AF 2 (R 2 ,s 2 ) ■ — 



h J v ' x ' ' \ h 
W{x 1 - ir 1 ,x 2 - ir 2 ,X 1 ,X 2 ,t;/3,X) 



(27) 
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Here we have defined full classical Liouville operator for the heavy mass subsystem as, 

P 8 8 
iLAXA = -L.^-+F R .--^—, (28) 
v ; M 8R l ' 8Pi' v ' 

for (i = 1,2), where the force F Ri = —dV{ji,Ri)jdRi. We have also introduced the quantity 

AFiiRi^i) = (v(r i ,R i )5(s i )- J dfV(n - f, RJ «>s(^)) . (29) 

If the potential is decomposed into light and heavy mass subsystem potentials, Vg(?*i) and Vh(Ri), respectively, and 
their interaction potential V c (ri,Ri) as V = Ve + Vh + V c , it is easy to demonstrate that AFi(Ri, Si) depends only on 
the interaction potential. 

The quantum-classical evolution equation (|27|l for W can be written formally and compactly as 

!L W ( t ) = ±KoW(t), (30) 

where the operator K is defined by comparison with Eq. I|27() . To simplify the notation, we have dropped the classical 
phase space arguments here and some of the following equations when confusion is unlikely to arise. 

This quantum-classical evolution equation for the spectral density is not yet in a convenient form for simulation 
since the kernels that appear in this equation are highly oscillatory functions arising from the fact that a Wigner 
representation of the quantum degrees of freedom has been used. In the next subsection we re-introduce the adiabatic 
basis and obtain a form of the quantum-classical evolution equation that can be solved by surface-hopping schemes. 



D. Quantum-classical evolution equation for W Q i aiQ 2 a 2 

The operators T and 1 1 can be used to convert Eq. (|30|l into an evolution equation for W, the matrix elements 
of W. Acting from the left with 1 1 on Eq. 130(1 and inserting unity in the form W = T o 1 1 o W = T o W, we 
obtain, 

= ^KoW(f), . (31) 

The transformed operator on right hand side of Eq. I|31|l can be calculated explicitly by a straightforward but lengthy 
calculation which is given in detail in Appendix [5] The result of this calculation is 
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xW ' M02 (Xi,X 2 ,t;/3,X) . (32) 

We see that the apparently formidable evolution operator 7C for the matrix elements of the spectral density, which 
depends on eight quantum indices, takes a simple form consisting of a difference of two quantum-classical Liouville 
operators, each acting separately on the classical phase space variables and quantum indices with labels 1 and 2. The 
quantum-classical Liouville operators are just those obtained in earlier derivations of the quantum-classical Liouville 
equation 



cias 
2], 



i^a^ai.p'^iiXt) — (iuj a > iai (Ri) + iL a ^ a .{X i )\8 a i,^ i 5 ai i3 i — J ' a >. au ^^(Xj) , (33) 

where uw(i?) = (E a (R) - E a >(R))/h and 

iL a > ai (X. t ) = g • A + 1 (f# (ifc) + (ifc)) • A , (34) 

is the classical Liouville operator involving the mean of the Hellmann-Feynman forces where Ffy = 
— {a; R\ 9Vv/ J^' R ) | a - ffy — — (a; R\ SH g^ R ^ \ a \ R)- Quantum transitions and bath momentum changes are described 
by 



Pi ( 1 9 \ 



Pi ( 1 9 \ 

--^r ' (l + ^S ai0i {Ri) ' ^ J , (35) 

where S ai p i = (E ai — Ep i )d ai p i (-£r ■ d^^)" 1 and d (U p ( —< af, R\ V^I/3,;; R > is the non-adiabatic coupling matrix 
element. 

The formal solution of Eq. (|32(1 is 

wa ' x a ia ' 2 a 2 ^ = f ex p(i/a)W(0) J . (36) 



'2 



Now that the quantum-classical evolution equation for the matrix elements of W and its formal solution have been 
determined, we can return to the calculation of the quantum-classical limit of the quantum time correlation function. 



IV. QUANTUM-CLASSICAL CORRELATION FUNCTION 



Equation l|32|l is one of the main results of this paper. Using it we can obtain a quantum-classical approximation 
to the quantum correlation function by replacing the full quantum evolution of the spectral density in Eq. Ill l[l with 
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its evolution in the quantum-classical limit given by Eq. ®, the solution of Eq. We have 

C AB (t;f3) = ~ f d\ J2 J fidX.B^iXMT' 2 ^) 

x(eM^t)W(X ll X 2 ,0;p,X)) . (37) 

Since the operator 7C is the sum of two operators, one acting only on functions of X\ and quantum indices with 
subscript 1, and the other on functions of X 2 and quantum indices with subscript 2, we may integrate by parts to 
have the operator act on the dynamical variables instead of W. We obtain, 

2 



where 



C AB (t,P) = £ /lI^^ 1 ^(X 1 ,|)^(X 2 ,-|)F^ /3l/3i,32 (X 1 ,X 2 ;/3) 



(38) 



a, Qi 



A^(X 2 ,-h = £ (e-^ x A R A^\x 2 ) . (39) 



a 2 a 2 

In writing Eq. |3*Hjl we defined 

W P ' 10M {X u X 2 ;p) = \ [ dXW^'^(Xi,X 2> 0;(3,X) ■ (40) 
P Jo 

Equation (|38|l shows that the correlation function at time t can be calculated by sampling Xi and A 2 from suitable 
weights determined by W 1 1 2 2 (Xi, A 2 ; /3) at time zero and propagating B^ 1 " 1 forward in time and A^ 2 -" 2 backward 
in time for an interval of length t/2. Note that while the time evolution in Eq. (|38|l is by quantum-classical dynamics, 
the initial condition for W^ 101 ^ 2 ^ 2 {X\, X 2 ; (3) is still an exact expression for the full equilibrium quantum mechanical 
spectral density. 



V. HIGH TEMPERATURE FORM OF W 



At t = W is given explicitly by 



W< a ^ a *(X u X 2 ,O;0,\) = 1 I ' dZ 1 dZ 2 e-^ p ^ +p ^ 

{2-Khy v h-z Q j 

x (a'^R^R, + ^\ e -^-^\R 2 - ^}|a 2 ; R 2 ) 

x (4;i? 2 |( J R 2 + ^| e -^|i? 1 -^)| ai ; J R 1 ) . (41) 
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It can be computed using path integral techniques but its evaluation is still a difficult problem. In order to illustrate 
its structure we consider its form in the high temperature limit. In this limit we may write 



Uh/2 M(R 12 - Z c f 



2Xh 2 



(42) 



where h = #- + V and we have introduced the variables Z c = [Z\ + Z 2 )/2, Z± 2 = Z\ — Z 2 , Rc = (Ri + -R2V2 and 
R12 = R\ — R 2 . Taking the desired matrix element of this expression and inserting complete sets of adiabatic states 
we obtain 

(a' 2 ;R 2 \(R 2 + ^\e7 xk \R x - ^)\av,Ri) = £ e~ XE ^-^ (a> 2 ; R 2 \a; R c - ^f)(a;R c - ^f\a 1 ;R 1 ) 

a 

= Y J ^ XEa{Rc) ^ l 2 ,R2W,R c ){a-R c \a 1 -R 1 ) + 0{Z 12 ). (43) 

Keeping only the zero order term in Z\ 2 , the integral over Z\ 2 in Eq. I|41|) gives a factor (An%) Vh 8(Pi — P 2 ). The other 
term in Eq. 141(1 . which arises from the combination of the gaussian on the right hand side of Eq. ((42(1 along with the 
analogous expression coming from the high temperature limit of + ^-\e~^~ x } H \R 2 — 4^), can be evaluated by 
performing the gaussian integral on Z c to obtain 

M \ Vh/2 j„ p 2A-ff D 2M p2 A(/3 — A) 2P? f ]\J \" h/2 A(3-A) 2P C 2 

\ e n2P^Ri2 e nW^ e —^r = _ f(R 12 ,P c )e — ^ , (44) 



2ttTi 2 I3) \2irh 2 f3 
where P c = (Pi + P 2 )/2. The function f(Ri 2 ,P c ) still contains quantum information since it is composed of a 
phase factor and a gaussian expressing quantum dispersion effects in the heavy mass coordinates. We can obtain a 
classical bath approximation if we represent f(R\ 2 ,P c ) in a multipole expansion and keep only the first order term, 
f(Ri2,P c ) « [f dR 12 f(R 12 ,P c )}S(R 12 ), we have 

/ £ 2 rt\ Vh l 2 p 2 

f(**>Pc) « Zr) e-**V»-fV5{R») (45) 



, 2M , 

Combining terms we obtain a high-temperature, classical-bath approximation to W: 



W< a ^(X u X 2 ,0;l3,X) « - * - e-P&e (P X)E < {Rl) e XE < {Rl h a , a J a , a , Rl 8{R 12 )5{P 12 ). (46) 

(2irn) Vh ZQ 1 2 



Thus the quantity W, defined in Eq. 1)40(1. is given in the high-temperature, classical-bath limit by 

W a ' iaia ' 2a2 (X 1 ,X 2 ;(3) = 



r( p i + P P(E a > (Ri)-E a > (Ri)) , 



x <5 Q;Q2 ^ a;fll ^(i?i 2 )(5(Pi2) . (47) 
Using similar manipulations, the high temperature limit of Zq is 
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If Eq. (|47|l is used in the correlation function formula, Eq. <|38[1 . the result maybe shown to correspond with the 
quantum-classical linear response theory form 29] to lowest order in ft. 

VI. CONCLUSION 

The expression for the quantum-classical limit of the quantum correlation function derived in this paper provides a 
route for the calculation of quantum transport properties in condensed phase systems. Difficult many-body quantum 
dynamics is replaced by quantum-classical evolution which can be carried out using surface- hopping schemes involving 
probabilistic sampling of quantum transitions, with associated momentum changes in the bath, and classical trajectory 
segments. The classical trajectory segments are accompanied by phase factors that account for quantum coherence 



when off-diagonal matrix elements appear. 21] The full equilibrium quantum structure of the entire system is retained. 
While the equilibrium calculation is still a difficult problem it is more tractable than the quantum dynamics needed 
to treat the many-body system using full quantum dynamics. For example, imaginary time Feynman path integral 
methods for computing equilibrium properties are far more tractable than their corresponding real time variants. Since 
quantum information about the entire system is retained in the equilibrium structure, the formula for the correlation 
function incorporates some aspects of nuclear bath quantum dispersion that is missing in other quantum-classical 
schemes. The importance of retaining the full quantum equilibrium structure has been noted in Ref. \&^ . 

The results also provide a framework for exploring and extending the statistical mechanics of quantum-classical 
systems. The correlation functions for transport properties that result from linear response theory in quantum-classical 
systems involve both quantum-classical evolution like that derived in this pap er, as well as the equilibrium quantum- 
classical density that is stationary under the quantum-classical evolution. 29J, |30j One may construct approximations 
to quantum transport properties by considering other approximate limiting forms of the equilibrium spectral density. 
We also note that to establish a complete comparison with quantum-classical linear response theory requires the 
retention of terms that were neglected in the calculations for W presented in Sec. [3 It should be fruitful to pursue 
extensions of such calculations to obtain other approximations for quantum transport properties. 
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Appendix A: DERIVATION OF QUANTUM-CLASSICAL EVOLUTION EQUATION FOR W 



The equation of motion for the spectral density (Eq. takes a similar form in scaled coordinates: 



-\ (*4 0) Vi, X'i) " l 4 0) '(4,^2)) W'ix'^x'^XiX^t'; (3', A') 

\ [ f[d S 'dS> (wUri.-i.^i^OW^)-^^,^^.^^)^) 



x W (x; - tt[ , X[ - n; , x' 2 — tt' 2 , X' 2 — II2, t'\ p', A') 



(Al) 



where the scaled free streaming Liouville operator and integral kernel are defined in Eqs. (|24(l and l|25|l . respectively. 
Inserting Eq. 1|26|1 into the expression for lo[ and retaining only terms up to linear order in in \i we find 



j[(r[, si, R[, S[) « J df'.dR'.V 1 ^ - f' 1 ,R' 1 ) sin fa ■ f[ + 2S[ ■ R[] 

. )n , R[ sin (2si • f[ + 2S[ ■ R[) + 0( t i 2 



W J 



(A2) 



We observe that 

<ifidRiV'(ri -fi,i?i)sin (2s; • f' 1 + 2S[ ■ R^j = j <lr\dR\V '{r\ - r\. l!\ ) sin (2.^ • r\ )v<»{2S' j ■ H] i 



(A3) 



+ 003(25; ■f' 1 )sm{2S[ ■ R[) 

using the trigonometric identity for the sine of a sum of arguments. Then, using the fact that J dR' 1 cos(2S' 1 R' 1 ) = 
7T Vh S(S' 1 ) we have 

df^dR^V (r[ ~ f£, i£) sin (2s[ ■ f[ + 2S[ ■ R[) = n Vb [ df[V' '(r'i - r[,R[) sm(2s[ ■ f^W]) , (A4) 



where we have used J dR[ sm(2S[ • R[) = 0. In a similar manner 



, ~, dV'(r[-f' 1 ,R' l ) 



df[dR[ 



dR[ 



R 



[ S mfa 1 -f' 1 +2Si-R! 1 ) = ~\ 



tt^ f J=f dV'(r[ - f[,R[) _ d5(S[) 



dr 1 



dR[ 



cos(2s 1 • r 1 ) 



dS[ ' 
(A5) 



where we have used the relations J dR\R\ cos(25^ • R\) = and / dR! x R! x sin(25*; • R[) = -(tt^ /2)dS(S' 1 )/dS[. Then 
to order 0(h), 



u'M^^Si) = — J d^V , (r[-f' 1 ,R[)siia(2s , 1 -f[)5(S , 1 ) 

^ y ^ — dR{ — cos(2si ' ri) ' ~ils[- 



(A6) 



Using this expression we may compute the integral on the right hand side of Eq. I|A1(I involving uj[ . The algebra 
for the uj' 2 term is similar. Given the expression l|A6() . the integral 



J ds\dS[J l (r[,s[,R l x ,S[)W\x[-^X[-U l 1 ,x' 2 ,X , 2 ,lf-,0,\) 



(A7) 
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has two contributions. The first is is 



-L J daldSiW'fa -n'^Xi- n' 1( 4, X 2 , f ; (3, A) J dfJ^Vi - r[,R[) sm(2s[ ■ r[)S(S[) = 

-L J d^W'to -^X[,x' 2 ,X' 2 ^,p,X) J df' l V'{r' 1 -f' 1 ,R , l ) S m{2s' l -f' 1 ) , (A8) 



while the second is 



J ds' 1 dS' 1 W'(x 1 -^,X{-n' 1 ,x' 2 ,X^t';(3,X) J dr[-^-V> (r[ - f[, R[) cos(2 s ' 1 • f[) • ^ 



7T 



(A9) 



Defining Ai* 1 ^, g , = — [V(i?^)(5(s' 1 ) — ^ J df^ cos(2si • r'^V'^ — f[, R[)\ and returning to unsealed coordinates 
we obtain Eq. I|27[) . the desired quantum-classical evolution equation for the full Wigner representation of W. The 
first term in the definition of AF' is compensated by the introduction of the full classical propagator for the heavy 
mass degrees of freedom. Written in this form, AF' also depends only on the interaction potential between the light 
and heavy mass particles. 



Appendix B: EQUATION IN THE PARTIAL WIGNER REPRESENTATION 



In this appendix we perform explicitly the calculations that, starting from Eq. (|27[) . lead to Eq. (|32[) . This calculation 
amounts to the evaluation of (T _1 o K o T) o W(t). The various terms in K, defined by Eq. I|27|l . are considered 
separately. 

Consider the calculation of (T _1 o iLi(X\) o T) o W which is composed of two terms. The force term is T _1 o 
Fri • d/dPiT o W = Ffjj • d/dP{W, since T and its inverse do not depend on Pi. The free streaming term 
(T _1 o -£l . JL- o T) o W requires additional calculations since T depends on R\. We have 



v M 



2 

x ^fe + y;i? 2 )0 Q2 (V 2 - y;i? 2 )^(7-2 + -^;Ri)4>% 2 {^ ~ y^ 2 ) 

j(r 1 + ^;ii 1 ) Zl Zl ^(n-f;^) 



dRi " PlVi 2 



OR! 



(Bl) 



2 '" 1J dR 1 

The last term where dW/dRi appears is simple to calculate and gives ■ ^-H/" 1 " 1 " 2 " 2 . To calculate the other 



two terms, we make the change of variables q± = r% 



-, Q2 = n 



f , 93 = r 2 



S 94 = r 2 



Integrating over 
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qs and q± and using J dq<fi a (q, R)<fi%(q, R) = 8 a p, we obtain 
Pi 



M 



P±P[ 



^2 / dqidq 2 <j}* a , i (q 2 ;Ri)(t) a A<li; R i) 



1 1) 6* 0i ( qi -R 1 ) + 6^( q2 ;R 1 ) 1 ' 



dRx 



dRi 



yy0' 1 Pia' 2 a 2 



(B2) 



Pi Pi 

where we have introduced the definition of the nonadiabatic coupling vector. 

Consider the calculation of (T _1 oiL^f^xi) oT) oW. In this case one integration over the momentum, arising from 
the definition of T _1 , gives J dp\{p\/m)e^ PllyZl ~ Z3 " 1 = {2^h) ui [ili/m)d8{z\ — z^j/dzs, while integration over p 2 gives 
(2irti) Ue 6(z2 — z^). One can then integrate by parts on z% and obtain 

2 3 

(T_1 °m ^"° T)oW = "S £ ffldr j f[dz k S(z 1 -z 3 )r a ' 1 {r 1 + ^;R 1 )ct >ai (r 1 -^-,R 1 ) 

1 PiPiPiP'J j=l k=l 



x0* ; (r 2 + ^;i? 2 )0 Q2 (r 2 - -|; ^2)^/3' (r 2 + -77; #2)^ (»"2 - -|;-R 2 ) 



y 2 



dz^dri 



siin + ^R^lin-T'^ 



(B3) 



Then use 



dz^dri 



,,{ n + ^-R 1 )^ i {r 1 -^;R 1 ) 



d 2 ^ l (r 1 + f;R 1 ) 



*3 



*3. d^^l-f^l) 



+ f 



p ^1 ( r i _ y 5 -Ri) - ^ (ri + y ; iZi)- _ a)2 



(B4) 



go back to the integral, perform the integration on Z3 and make the change of variables previously introduced. The 
integration over q% and 94 can be performed using the completeness of the adiabatic basis and one gets 



(T- 



-1 Q Pi Jl_ 
to dr\ 

i 



oT)oW= l - y2(a' l] R 1 \^-\P' 1 ;R 1 }W^ a ^ 
n 2m 



Pi 



^(/3 1 ;i? 1 ||-|a 1 ;i? 1 }^' 3l ^ Q2 , 
Pi m 



(B5) 



where we have used the identity (a, R\p 2 /2m\P, R) = -(h 2 /2m) J dq6* a {q 1 R)d6 (q^R)/dq 2 . 

We consider now the transformation of the potential term in K o W equal to c\ J dsids 2 8(s 2 ) J dfV(r\ — 
f ) sin (2sp) W {xi - m,x 2 - TT 2 ,X 1 ,X 2 ,t;PX), where c x = 2h~ 1 (7rfi)-"«; we may immediately perform the trivial 
integrations on pi, p 2 , z 3 and z 4 . We also make the change variables below Eq. 1|B1|) so that the integrations on q 3 
and (?4 are also trivial. One is left with the integral 



P[Pi 



Ci / dqidq 2 / dsidrV( 



r) sin 



«i(?2-gi) 



x6* a{ (q 2 ;R l )4> ai {qi;Ri)<t>P[ faRiWfh foil Ri)W^ a > a > 



(B6) 
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Using the fact that J ds 1 e^ Sl ' < - q2 - qi) sin (2si ■ f/%) = {2ith) Vl [£(g 2 - <7i + 2f) - <5(q 2 - qi - 2f)] /2i, substituting into 
Eq. I|B6|) and making the change of variable a = 2f, the delta functions can be integrated out and one obtains, 



-y 



ih 



dqi<f> ai (qi; Ri)(f>* 01 (qi;Ri) 



dq2 C; (92 ; Ri ) V(q 2 ) (frft (q 2 ; Ri ) 



■4 I d qi ( l ) p 1 (q^ R i) v (qi) < l > c ll (qi;Ri) f (92; #1)^(52; -Ri) 

* ftft 1,7 J LJ ' 



= 4 V{ai; J Ri|t>|^i; J Ri)W 3 ^ ia = a2 - — V(/3i; i?i|F|ai; i?i)W^ /3ia ^ . (B7) 
i7i z — ' in z — ' 

/3; ft 

The last term that must be worked out explicitly arises from the transformation of the quantum-classical term 
J dsids2S(s 2 )AFi(Ri,si) ■ -^W{xi — Wi,x 2 — ir 2 , X\, X 2 , t; A, (5). Recalling the expression for Ai*\ in Eq. (|2^|) one 
sees that there are two contributions to transform. The transformation of the term involving 9 g^ 1 ' 1 <5(si) is the same 
as that for the force term in iLi(Xi) and yields Fji t ■ dW aiaia2a2 / dP\. The integral term in Aiq(i?i,si) can be 
computed by integrating over pi, p 2 , Z3 and 24, performing the change of variables below Eq. (|Blf) and integrating 
over (73 and 94, to obtain 



y l %\u e [ i(q2i R i)0a 1 (qi;Ri)<t>0[(q2iRi)<t>X( , ii. : > R i) 

a/ a. y Kl% ! J 



0[P 

x / dsi I df 



V( — r, Ri) ] cos 



dRi 



2si • f 



d 



(B8) 



Then one can use the integral / ds 1 cos (2s x • r/K) e ^ Sl ' (q2 ^ qi '> = {2%%) Vt [8{q 2 - q 1 + 2f) + S(q 2 - q 1 - 2f)] /2, make 
the change of variable a = 2f, and integrate out the delta functions to find 



ft 



(B9) 



Analogous terms (but with opposite sign) are obtained when considering the transformation of the terms depending 
on x 2 and X 2 in Eq. f2Tjl . 

Combining all terms and using the relations F^{R) = -(a; R\V R V(R)\f3; R) — F$ + (E a - E )d a p, 



Pi 



^(ai;^! (t^ + v) \ft;R 1 )W*°'^'*> + + \a i; R^W^' 2 * 2 = 



2m 



{E ai {R{) - E a ,{R x ))W a >^ a * = -iu a{ai (Ri)W< aia !> a > 



(BIO) 



and introducing the definition S ai /3 i — (E ai — Ep^da^p^^ • d^fa) 1 , we find Eq. I|32|) 
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